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Ch Abstract 

^~5 Interest in continuous-time processes has increased rapidly in recent years, largely because 

\^ of the high-frequency data available in many applications, particularly in finance and turbu- 

lence. We develop a method for estimating the kernel function of a stationary Levy-driven 

L^ continuous-time moving average (CMA) process Y based on observations of the discrete-time 

r^ process Y^ obtained by sampHng F on a discrete grid. We define a non-parametric estima- 

tor of g^ , based on the Wold representation of Y^, and show that for every CARMA(p, g) 
process g^ converges pointwise to the kernel g of the process F as A 4- 0. Since the esti- 
mator is non-parametric, it is well-defined for all CMA processes with absolutely continuous 
spectral distribution function. We illustrate its performance by applying it to simulated high- 
frequency data from a CARMA process and a process with non-rational spectral density, and 
C^ then to the Brookhaven turbulent wind speed data. In the course of the analysis we extend 

^ the results of Brockwell et al. (20121 to a wide class of CMA processes with strictly positive 






00 

Kf^ spectral densities, and we derive higher order asymptotic approximations, which are valid 

■^ for CARMA processes of all orders. 

^, 
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1 Introduction 



We are concerned in this paper with causal continuous-time moving averages of the form 

/oo 
g{t-s)dLs, te^, (1.1) 

-oo 

where {Lt]t&. is a Levy process with ELi = and EL\ = a^ < oo. The kernel function g is 
assumed to be zero on (— cx),0] (for causality) and square integrable, but it may be singular at 



the origin as in Example 4.8 The process Y defined by (1.1) is then a zero-mean strictly and 



weakly stationary process. For the estimation of g discussed in Sections 5 and 6 we make the 
additional assumption that L has been standardized so that EL^ = 1, since otherwise g and cr^ 
are confounded. 

The integral in (1.1) is understood in the L^-sense and, since we use only second-order 



properties in our analysis, the results apply more generally to processes defined by (1.1) with 

and ELl = a"^ as in |Doob| (|1990"1) 



L a stationary orthogonal increment process with ELi = u ana J^i^j = a" as m 
Ch. IX. It is important to note however that, when L is a given Levy process, Y is completely 
characterized by g^ while the spectral density of Y characterizes only the second-order properties. 
Throughout this paper, stationarity will always mean weak stationarity. 

Examples of CMA processes are the Ornstein-Uhlenbeck process, with g{t) = e'*'*l(o^oo)) where 
A < 0, and the more general continuous-time autoregressive moving average (CARMA) processes 



studied by Doob (1944) for Gaussian L. State-space representations of these processes were 



exploited by Jones ( 1981[ ) and Jones and Ackerson (1990) for dealing with missing values in time 
series, and by Brockwell (2001) for the study of Levy-driven CARMA processes. Long-memory 



versions have been developed by Brockwell and Marquardt (2005) and Marquardt (2006). CMA 



processes constitute a very large class of continuous-time stationary processes (cf. Doob 



(1990) 



Theorem 5.3, Ch. XII and Yaglom (2005), Section 26.2). A gamma kernel function g corresponds 



to a process with Whittle-Matern autocovariance function and spectral density of similar form to 



the von Karman and Kaimal spectral densities, which have been widely employed in turbulence 



modelling. The extremal properties of CMA processes have been studied in Fasen (2009). 



The work of the present paper was originally motivated by a study of the Brookhaven turbu- 



lence data (see Ferrazzano (2010) for a detailed description). The data consists of twenty million 
values, sampled at 5000Hz (i.e. 5000 values recorded per second over a time interval of 4000 
seconds). In order to understand the relationship between such high-frequency data and the un- 
derlying continuous-time process Brockwell et al. ( 2012[ ) studied the asymptotic behaviour of the 
sampled process Y = (l^A)nez as A ^ when 1" is a CARMA(p, q) process with p — g < 3. 

In Section 2 of this paper we derive higher-order asymptotics for Y which apply (unlike 
those of Brockwell et al. (2012 ) to all CARMA(p, q) processes and we use the results of Section 2 



in Section 3 to establish the pointwise convergence of a family of functions g( , defined in terms 
of the Wold representation of Y , to 5 as A J, 0. In Sections 5 and 6, we use non-parametric 
estimators of g to estimate g based on observations of Y with A small and illustrate the 
performance of the estimators using both simulated data and the Brookhaven series described 
in the preceding paragraph. The outcome of a detailed statistical analysis for turbulence data is 



presented in Ferrazzano and Kliippelberg (2012). 



In Section 4 we extend the asymptotic results of Brockwell et al. (2012 ) in a different direction 



by deriving analogous results for a broader class of CMA processes with strictly positive spectral 
density. 



We use the following notation throughout: ?R.{z) denotes the real part of the complex num- 
ber z; B denotes the backward shift operator, BY^ := Y^^ for n G Z; a(A) ~ b{A) means 
lim^io o(A)/6(A) = 1. As g can have a singularity in 0, the spectral densities of Y or Y may 
have a singularity in as well, and we may have to restrict the range of frequencies for their 
spectral densities to Qc '■= ^\{0} and Q^ := [— 7r,7r]\{0}. 

2 Asymptotic behaviour of Y^ as A ^ 



The CMA process Y defined by (1.1) has autocovariance function 



^Y{h) = cr^ / g{x)g{x + h)dx, /i G 



and spectral density 



where 



/yH 



1 

2^ 



-iujh 



-f{h)dh 



a 

2^ 



|5(w)P, w G ilc 



(2.1) 



5(w) : = 



'^g{x)dx. 



The spectral density of the sampled process, Y := iXn/\)n&. is (Bloomfield (2000), p. 196, 
Eq. 9.17) 

'uj + 2A;7r^ 



k=—oo 



A 



io G VLa. 



(2.2) 



For the causal, finite variance Levy-driven CARMA(|?, q) process, autoregressive polynomial 
ai^z) and moving average polynomial 6(2), the spectral density is 

2 



/y(^) 



a 

2^ 



h{iiS) 



ayiLo) 



-TT < a; < vr, 



(2.3) 



where 0(2:) = z^ ^a\z^ ^ -|- • • • -l-Op, 6(2;) = 6o + ^i-2 + - • • ^^qZ^ ^, p > q, bg = 1, and the zeros of 
a{z) all have strictly negative real parts. Without loss of generality we can also assume that a{z) 



and b{z) have no common zeros (see Brockwell and Lindner (2009), Theorem 4.1). The kernel is 

Hz) 

-^e uz ±(o^oo){i) = / ^ixt3S2=A \ K 



^(^) = i/^^"^^^(°-)(*)=^^^^^ 



a(z 



l(0,oo)(*), 



(2.4) 



where the integration is anticlockwise around any simple closed curve p in the interior of the left 
half of the complex plane, encircling the distinct zeroes A of a{z), and ^eSz=\{f{z)) denotes the 



residue of the function / at A. For such processes it was shown in Brockwell et al. (2012 ), Section 
2, that the spectral density /a of the sampled process Y is 

, , , -a"^ f b(z)b(-z) sinh(A2) , , , 

^^^") = « y,a(z)a(-z)cosh(A.)-cos(.)^"' "^ ^ ^ ^ -' (^.5) 



where the integral, as in (2.4), is anticlockwise around any simple closed contour p in the interior 



of the left half of the complex plane, enclosing the zeroes of a{z). It is known that the sampled 
process Y satisfies the ARMA equations. 



^a{B)Y^ = eA{B)Zt nGZ, {Z^}nez^WNiO,ai) 



(2.6) 



where B is the backward shift operator, 9/^{z) is a polynomial of degree less than p, [Z^)nez is 
an uncorrelated sequence of zero-mean random variables with variance, which we denote by u^, 

p 

Mz) = ll{l-e^^^z), zee, 



and Ai,...,Ap are the zeroes of the polynomial a{z). Since the polynomial 0a (-z) is known 
precisely for any given CARMA process, the second-order properties of the sampled process 
Y for small A can be determined by studying the properties of the moving average term, 
Xn '■= 9i^{B)Z^ in (2.6), as A ^ 0. Denoting by fuA the spectral density of X, we find from 



(2.6) that 



/maH = 2Pe-'^i^/A(^) n(cosh(A,A) - cos(l^)). 



-vr < (J < TT. 



(2.7) 



Brockwell et al. (2012) determined the leading terms in the expansions of /a and fuA 



in powers of A. These terms determine the local second-order behaviour of the corresponding 
processes. In Section 4 we extend these results to a more general class of CMA processes. 

In the following section we introduce a small-A approximation g to the kernel g oiY based 
on the Wold representation of the sampled process Y . In order to show the convergence of g 
to 5( as A ^ for CARMA(p, g) processes, we need to consider higher order expansions of the 



spectral densities /a and Jma than were considered in Brockwell et al. (2012). We conclude this 
section by deriving the required expansions. 



From (2.5) it follows at once that the spectral density fA.i^) is — cr^/(27r) times the sum of 
the residues at the singularities of the integrand in the left half-plane, or more simply o"^/(47r) 
times the residue of the integrand at oo, which is much simpler to calculate. Thus, 



fAico) 



Att' 



-Res, 



biz)bi-z) 



sinh(A2;) 



a{z)a{—z) cosh(A2;) — cos(a;) 



-TT < CO < TT. 



The spectral density can also be expressed as a power series. 



2 

/a(^) = ^ ^a'A2j+V,c,(a;), -Tr<u<7T, 



{21 



j=0 



where Cfc(w) is the coefficient of z ~^^ in 



^Cfc(a;); 



,2fc+l 



sinhz 



A;=0 



cosh z — cos u) 



-vr < cj < vr. 



and 



r, := ReSz 



gy+i Kz)b{-z) 



a(z)a[—z) 
i.e. the coefficient of z"^^ in the power series expansion. 



Y^r.z'^ = i-z') 



2^2^ 



2\p-q-l i[i=l\^ f^jZ 



j=0 



nLi(i-A?-^ 



(2.9) 



where a{z) = YVi=i{z — Xi) and b{z) = Yl1^-^^{z — m). The power series (2.8) is the required 



expansion for /a- The expansion for Jma is obtained from (2.7) and (2.8) as 

i=l 



hlAi^) 



Att 



n(i 



cosw 



+f;^)Ew-)A--\ -.<^<^. 



t (2J) 



j=i - - - fc=0 

This can be simphfied by reexpressing it in terms of x := 1 — cosw. Thus 



^^-^^"^= 4vr n(^ + E (2,), 



(A.A)2^- ^ 



y^rfcafc(x)A 



2A:+1 



(2.10) 



1=1 j= 

where afc(x) is the coefficient of z^^^'^ in the expansion, 

/c=0 



fc=0 



cosh z — 1 + X 



In particular 00(2;) = l/x,ai(x) = (x — 3)/(3!x^) and a2(a;) = {x^ — 15x + 30)/(5!x'^). More 
generally, afc(x) has the form. 



afc(x) 



(2A; + l)!x'=+i JLi 



n(^~^'='*)' 



where 



i=l 



(2.11) 



(2.12) 



and the product, when /c = 0, is defined to be 1. Since ap-q-i{x) plays a particularly important 
role in what follows, we shall denote its zeroes more simply as 



si • — sp— g— l,i) * — -I-) • • • )P 



1. 



From (2.10), with the aid of (2.9) and (2.11), we can now derive the required higher-order 



approximation to fMA{^)- Observe first that the expression on the right of (2.10), in spite of its 



forbidding appearance, is in fact a polynomial in x of degree less than p. We therefore collect 



together the coefficients of x^ ^,xP ^, . . . ,x'^. This gives (using the identity (2.12) and defining 
y := A^) the asymptotic expression as A 4, 0, 



fMA{(^) 



2Po-2g-aiA^2(p-g)-l 

47r 



x^r. 



p_q_lQp_y_l(x) +0(1) +^PjX^ ^y^ 



(2.13) 



with 



Pj 



(_2)-(p-.-i+:.-) [r,^,_,+, - rp_,_2+, f] A,^] + 0(1) 



i=l 



2-(p-.-i+.)j;^^...;,2^+o(l), 



where the second line follows from ( |2.9[ ) and the sum on the second line is over all subsets of size 
j of the q zeroes of the polynomal b{z). 



Finally, replacing Vp-g-i in (2.13) by (—1)'*'^'^^^, substituting for ap-q-i{x) from (2.11) and 
using the continuity of the zeroes of a polynomial as functions of its coefficients, we can rewrite 
(2.13) (recalling that x := 1 — coslo and £^i,i = 1, . . . ,p — q — 1 are the zeroes of ap_q_i(x)) as 

9 






[2{p - q) - IYAtt 
Observe now that we can write 

4 



,2A2 



«=1 



1 



k=l 



//2A2 

:r+ ^(1 + 0(1)) 



X + ^^:^(1 + oil)) = —(1 - ae-^")(l - Cjen, 



-IT < UJ < TT, 



where 



a = l±^fcA + o(A), 
and the sign is chosen so that limA4.o iCfcl < 1- Similarly we can write 



(2.14) 



(2.15) 



x-e,(l + o(l)) 



M^i) 



[l + v{^^)e-n{^ + v{^^)en, 



where 



rji^,) =Ci-l± Vte -1)2-1 + 0(1) 



(2.16) 



and the sign is chosen so that liniAj^o I^(Cj)I < 1- If the zero S,i of ap-q-i{x) is such that both 
choices of sign cause the limit to be 1, then either choice will do provided the same choice is 
made for 77(Cj), where ^j denotes the complex conjugate of ^j . 

These factorizations allow us to give the following asymptotic representation of the moving 
average process X„ = 9i^{B)Z^ appearing in (2.6). 

Theorem 2.1. The moving average process {Xn}neZ with spectral density Jma has the asymp- 
totic representation, as A ^0, 



p-i-g 



Xn= II {l + v{Ci)B)l[{l-CkB)Z^, {Z^}„ez~WN(0,ai) 



i=l 



fc=l 



where 



4 



^2(p-g)-lg-aiA^2 



[2{p-q)-lv■mr'v{^^)nuck 

with Qk CLnd r]{^i) as in (2.15) and ^2.16^. 



Proof. The result follows at once from (2.14), (2.15) and (2.16). 



(2.17) 



(2.18) 



D 



Remark 2.2. (i) The parameters r/(^i) and (^k may be complex but the moving average operator 
will have real coefficients because of the existence of corresponding complex conjugate parameters 
in the product. 



(ii) The representation in Theorem 2.1 is a substantial generalization of the one in Corollary 2 of 



Brockwell et al. (2012 ), since it is not only of higher-order in A, but it applies to all CARMA(p, i 



processes, not only to those with p — q < 3. 



a 



3 The Wold approximation to the CARMA(j9, g) kernel 

In this section we introduce an approximation g to the kernel g of the process Y defined in 



(1.1 1, which depends only on the Wold representation, 






(3.1) 



of the sampled process Y . The approximation is 



i=o ^ ^ 



(3.2) 



Using Theorem 



2.1 



we shall show that, for all CARMA(p, q) processes, as A ^ 0, ^f converges 
pointwise to og^ or to gf if L is standardized so that EIj\ = 1. We first illustrate the convergence 
in the simple case when y is a CARMA(1,0) or stationary Ornstein-Uhlenbeck process. 



Example 3.1. [The CARMA(1,0) process] This a special case of (1.1) with kernel 



9\x) 



e^^l 



(0,oo) 



ix) where A < 0. 



The sampled process 1" is the discrete-time AR(1) process satisfying 

where Z = {^„ } ^^ is the independent and identically distributed sequence defined by 

r-nA 



Zt 



gA{nA-«)^^^^ n G Z. 



(n-l)A 



In this case it is easy to write down the coefficients 'i/'^ and the white noise variance a\ in the 
Wold representation of Y . From well-known properties of discrete-time AR(1) processes, they 
0, 1,2,.. ., and a\ = %^{e'^ — !)• Substituting these values in the definition 



areV'f = e^^^,j 



(3.2) we find that 



g^'ix) 



which converges pointwise to ag as A ^ 0. 



,2AA 



pjAA 



2AA 



e-' ^jA,{j+i)A){x), 



a 



The approximation (3.2 ) is well defined for all processes ( 1.1 ) and there are standard methods 



for estimating the coefficients and white noise variance appearing in the definition from obser- 
vations of Y . Example 3.1 shows that g converges pointwise to ag for CAR(l) processes. 
Our aim now is to establish this convergence for all CARMA(p, q) processes. We give the proof 
under the assumption that the zeroes Ai, . . . , Ap of the autoregressive polynomial a{z) all have 
multiplicity one. Multiple roots can be handled by supposing them to be separated and letting 
the separation(s) converge to zero. 



The kernel (2.4) of a causal CARMA(p,q) process Y whose autoregressive roots each have 



multiplicity one reduces (see e.g. [Brockwell and Lindner (2009)) to 



5(^) = 2^ ;w 



-e^'^l(o,oo)(a:) 



(3.3) 



where a{z) = J|f^]^(2; — Aj) and b{z) = Y\1=iiz — m) are the autoregressive and moving average 
polynomials respectively and a' denotes the derivative of the function a. We now establish the 
convergence, as A ^ 0, of g as defined in (3.2) to ag. Theorem 2.1 is used to determine the 



parameters of the Wold representation appearing in the definition of g . 



Theorem 3.2. IfY is the CARMA{p,q) process with kernel (3.3), 

(i) the Wold coefficients and white noise variance of the sampled process Y are 



^ 



p 

E 

r=l 



nr=i"^(i+r/te)e-^'-^)nLi(i-ae- 



ArA\ 



JKA 



and 



-I 






[2(p-'z)-l]!^£^S(c^)^LlC/ 

with (k CL^d rj{^i) as in (2.15) and \2.16) and 



(3.4) 



(3.5) 



(a) the approximation g defined by (3.2) with ip^ and a\ as in (3.4) and (3.5) converges 
pointwise to ag with g as in (3.3). 



Proof, (i) The expression for o"^ was found already as part of Theorem 
is the coefficient of z^ in the power series expansion, 



2.1 



The coefficient ipf' 



j=0 



n£i"'^(i+^te)^)nLi(i-Cfc^) 



ui.=ia 



gAmA--") 



which can be seen, by partial fraction expansion, to be equal to (3.4). 



3.2[ ), substituting for (^ from (2.15), letting A 4- and using the identities 

a {Xr 



ii) The convergence of g follows by substituting for tj;^ and a"^ from (3.4) and (3.5) into 

Am) 



and 



p-q-l 

n 



:i+^te))^ 






n I = [2(p -(?)-!]!, 



the last equality following from (|2.12). 



D 

Remark 3.3. Although we have established the convergence of g only for CARMA processes, 
the non-parametric nature of g strongly suggests that the result is true for all processes defined 
as in (1.1). In practice we have found that estimation of ag by estimation of g with A small 
works extremely well for simulated processes with non-rational spectral densities also. □ 



4 Asymptotics for a class of sampled CMA processes as A J. 



Brockwell et al.] ( |2012[ ) derived first-order asymptotic expressions, as A ^ 0, for the spectral 
density /a when K is a CARMA(p, q) process with p — q < 3. Although, as pointed out in 
Section 2, these asymptotic expressions are not sufficiently precise to establish the convergence 



of g to g, they do reveal the local second-order behaviour of the process Y. For example, if Y is 
a CARMA{p,p — 1) process driven by a Levy process L with Var(Li) = ct^ then equations (15) 
and (19) of |Brockwell et al.| ( poT2| ) give, as A | 0, 



/a(^) 



47r(l 



cosw 



-7T < UJ < TT, 



showing that the spectral density of the normalized dijferenced sequence {{YnA—Y(n-i)A)/yA.}n£z 
converges to that of white noise with variance o"'^ as A ^ 0. In other words, for any fixed positive 
integer k, the sequence of observations 1^,a/vA, n = l,...,/c, from a second-order point of 
view, behaves as A ^ like a sequence of observations of integrated white noise with white-noise 
variance a'^. 

In this section we derive analogous asymptotic approximations for the spectral densities of 
more general CMA processes and the implications for their local second-order behaviour. Since 
we allow in this section for spectral densities with a singularity at zero we recall the definition 
of the spectral domains, 

Qd := [-7r,7r]\{0} and Q^ ■= (-oo,oo)\{0}. 

We require the CMA processes to have spectral density satisfying a weak regularity condition at 
infinity. To formulate this condition we first need a definition. 

Definition 4.1 (Regularly varying function (cf. Bingham et al. ( 1987| )). Let f be a positive, 
measurable function defined on (0, oo). // there exists p S M such that 



/(Ax) 



lim 



XP for all A > 0, 



holds, f is called a regularly varying function of index /) at cxd. The convergence is then auto- 
matically locally uniform in A. We shall denote this class of functions by TZp[oo). Furthermore 
we shall say that /(•) G lZp{Q-\-) if and only ?//(!/•) S TZ-p{oo). 



The characterization theorem for regularly varying functions (Theorem 1.4.1. in Bingham 
et al. (1987)) tells us that / S lZp{co) if and only if f{x) = xPL{x), where L £ TZo{oo). 



Theorem 4.2. Let Y be the CMA process ( |1.1[ ) with strictly positive spectral density fy such 
that fy G TZ^a{oo), where a > 1, i.e., for L £ TZo{(X)), 



\U}\ 



'L(|a;|), u £^c 



(4.1) 



/y(^) 

Then the following assertions hold. 

(a) The spectral density of the sampled process Y has for A ^ the asymptotic representation 



/a (a;) ~ L(A-i)A"-i [|a;|- + (2^)-"C (a, 1 " ^) + (2vr)-"C («, 1 + 



CO 

2^ 



CO £ Qd 



(4.2) 



where (^(s, r) is the Hurwitz zeta function, defined as 



Cis,r):=Y, 



fc=0 



{r + ky 



SR(s) > 1, r/ 0,-1, -2, 



(h) The right hand side of (4.2) is not integrahle for any A > 0. However, the corresponding 



asymptotic spectral density of the differenced sequence (1 — B)"'''^Y is integrahle for each fixed 
A > and the spectral density of 



[I - Br/^ 



L(A-i)i/2A("-i)/2 



Y' 



(4.3) 



converges as A ^ to that of a short-memory stationary process, i.e. a stationary process with 
spectral density hounded in a neighhourhood of the origin. 



(c) The variance of the innovations {Z^}n£Z ^'^ the Wold representation (3.1) ofY satisfies 

ai ~ 27rC„L (A-i) A°-S A i 0, 
where 



Ca = exp <; ^ / log 



IL^I " + (2^)-C(a,l-^l + 



(27r)-"C (a, 1 + ;^ ) 



2-kJ 



duj 



(4.4) 



Proof, (a) Since /y is positive, Eq. (2.2) can be rewritten as 



/.M = A-VK(A-)f;M^±M^. .e«.. 



fc=— oo 



/y(A- 



(4.5) 



Each of the summands converges by regular variation to |a; + 2TTk\ ". It remains to show that 
we can interchange the infinite sum with this limit. Invoking the Potter bounds (Theorem 1.5.6 



(iii) of Bingham et al. (1987)), for every e > there exists a A^, such that for all A < A^ and 

/y(|L^ + 27rA;|A-i: 



\2TTk + uj\ > 



(l-e)|27rA: + a;|-"~^ < 



/y(A^ 



< (1 + e)\2TTk + uj\ 



(4.6) 



We take e > such that a — e > 1. Then, using (4.6), we can bound (4.5) as follows: 
(1-e) ^^^^ ^ Y, |2^fc + t^|-"-^</A(c^)<(l + e) -^^^^ ^ Yl |27^A: + a;|-"+^ uj e ^a. 



fc=— oo 



fc=— oo 



(4.7) 



Since e can be chosen arbitrarily small, we conclude that as A J, 

— 1 \ oo 

a" 



fAioj) - ^"^^^ ^ E |w + 2A:vr|-°, c^ G O^. 



fc=— oo 



We can rewrite the sum above as 



y Iw + 2A;7rr" = (27r)-" V — + A; 

fc=— oo fc=— oo 



\LU\ 



- + (2vr)-"X^ 



k=0 



k + 1 + {k + l-\ 

2ttJ \ 2irJ 



(4.8) 
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From this and the definition of (^ we obtain (4.2). 



(b) We first note that the Hurwitz zeta function Ci~Ci^ Huj/Itt) is bounded and strictly positive 
for all u € 0^, therefore, its integral over [— vr, tt] is positive and finite. On the other hand, 
since a > 1, the term uj~°^ is not integrable over [— 7r,7r]. However, the differenced sequence 
(1 — i3)°/^y^, has spectral density 



h^ioj) = 2"/2(i _ cos6<j)"/Va(w), uj £ ^d. 



As A 4, we can write, for uj G fid, by (4.2) 



(4.9) 



h^iu) 



2°/2(i _cosa;)°/2L(A-i)A"-^ x 

CO \ 



\UJ\ 



+ (27r)-°C ( "' 1 - f^ j + (2vr)'"C [a, 1 + ^^ 



UJ 



The right hand side is integrable over [— vr, vr] and bounded in a neighbourhood of the origin, 
since 2"'^(1 — cosa;)°'^Ct;~" — )■ 1 as cj — )• 0. Thus we conclude that the spectral density of the 



rescaled differenced sequence (4.3) converges to that of a short-memory stationary process. 



(c) It is easy to check that the sampled CMA process has a Wold representation of the form 
(3.1) and that its one-step prediction mean-squared error based on the infinite past is a\. Kol- 



mogorov's formula (see, e.g.. Theorem 5.8.1 of Brockwell and Davis (1991)) states that the one- 
step prediction mean-squared error for a discrete-time stationary process with spectral density 
/is 



r = 2vr exp 



1 

2tt 



^og f{uj)du} 



(4.10) 



Applying it to the differenced process we find that its one-step prediction mean-squared error is 

27rexp< — / log/i (c<j)(ia; 
2vr ,1^ 



ol. 



= 2vrexp<-— / log(2 — 2cosa;)dLi; > x exp < — / log f^{uj)doj 
I4vr7__^ J I2vr7_^ 

Hence the differenced sequence has the same one-step prediction mean-squared error as Y itself. 



Since from (4.7), as A | 0, 

log/A(u;)-log(L(A-i)A"-i) 



log 



00 



2'Kk + uj\ 



pointwise on Vtd^ and since the left side is dominated by an integrable function on ^4, we conclude 
from the dominated convergence theorem that, as A J, 0, 



L(A-i)A°-i''''Pl2vr 



\ogf^{u:)duj > — )• exp 



2vr 



log 



^\2'Kk + u:\~'^duj 



which, with (4.8) and (4.10), shows that as A ^ 0, 

di ~ 27rC7„L (A-i) A"-^ 



(4.11) 



□ 
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Remark 4.3. (i) Theorem |4.2[ b) means that, from a second-order point of view, a sample 
{Y^,n = 1, . . . , fc} with k fixed and A small resembles a sample from an (a/2)-times integrated 
short-memory stationary sequence. If in (b) we replace (1 — i?)°'^ by {l — Bp where 7 > (a — 1)/2, 
then the conclusion holds for the overdifferenced process. If, for example, we difference at order 
7 = [(a -|- l)/2j (the smallest integer greater than (a — l)/2) we get a stationary process. In 
particular, if 1 < a < 3, then [{a + l)/2j 

(1 



1 and, by (4.2) and (4.9), the differenced sequence 



L{A- 



LJ (^ Qq 



B)Y has the asymptotic spectral density, as A J, 0, 

\"-i2(l-cos..)[|a;|- + (2^)-"C(«,l-^) , .-. , ^^^, ^ , 2vr. 

This is the spectral density of the increment process of a self-similar process with self-similarity 
parameter H = (a — l)/2 (see Beran (1992), eq. (2)). Moreover, for a generic a > 1 the asymp- 



+ (27r)-"C a,l + 



totic autocorrelation function of the filtered sequence has unbounded support. The only notable 
exception is when a is even, where the asymptotic autocorrelation sequence is the one of a 



moving-average process with order a/2, as in Brockwell et al. (2012) or in Example 4.7 



ii) The constant C^ of (4.4) is shown as a function of a in Figure 1. The values, when a is an 



even positive integer, can be derived from (3.5) since CARMA processes constitute a subclass 



of the processes covered by the theorem (see Example 4.7). It is clear from (4.4) that Ca is 
exponentially bounded as a — ?• 00. □ 



Coefficient C^ 



Coefficient d (log scale) 



O 



5 10 15 20 25 30 



Figure 1: The constant Ca, as a function of the index of regular variation a, is shown on the 
left using a linear scale and on the right using a logarithmic scale. From Corollary 3.4 (a) of 



Brockwell et al. (2012) we know that C2 = 1. The horizontal line indicates the value 1. 



Corollary 4.4. Let Y be a CMA process satisfying the assumptions of Theorem 4-2 with 1 < 
a<2p+l. Then for A 10, 

E[((l - BfY^f] ~ 2PSp,aL{A-')A"-\ 



where 



S, 



p,a 



fl-coswF 



\ur + (2^)-"C ( «' 1 - ^) + (2^)""C («' 1 + 1^ 



dtj. 
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Proof. By stationarity we have E[(l — ByY^] = and, hence E[((l — B)PY^)'^] is the variance, 
of ((1 — B^Yj^) which can be calculated as the integral of its spectral density. Thus 



E[((l - B)PY^y] = 2P (1 - cosuj)PUioj)du; 



Using the inequalities (4.7) and Lebesgue's dominated convergence theorem, we find that as 
A 4 0, 



L(A- 



——3^/ {1- COS uj)PfA{oj)duj^ {1 - cos u})P "^ \2kTr + uj\-'^du} 



fe= — oo 



which, with the previous equation and (4.8), gives the result. 



n 



The kernel of the CMA process (1.1) and its spectral density are linked by formula (2.1) 



Moreover, it has long been known that local properties of a function imply global properties of 



An Abelian theorem of Cline 



its Fourier transform (see e.g. Titchmarsh (1948), Theorems 85 and 86) 



(1991 ) allows us to show, under the conditions of the following 



proposition, that CMA processes with regularly varying kernels at the origin have regularly 
varying spectral densities at infinity. 

Proposition 4.5. Let Y he a CMA process with kernel g G TZu-i{0+) for v > 1/2. Assume that 
the derivatives in satisfy the assumptions 

(Al) 5(M)(o+)/0; 

(A2) 5(L''-iJ) E 7l«(0+) for a G [0, 1) (with g^-^^ := jlg{s)ds); 

(A3) For some xq > 0, 



[u) := sup sup 

x<xo 0<'w<v<l 



is bounded and integrable on [l,oo) 
Then 



gil''-^\){{u + V + w)x) - c/(L!--iJ)((u + v)x) - 5(L'^-iJ)((n + w)x) + g^^^-^i^ux) 



g(l'^-'^\){x) 



/y(| • I) G 7^_2.(oo). 



Proof. Under conditions (A1)-(A3) we can apply Theorem 2 of Cline (1991), which yields 



F{g){\u\) ~ T{u + l)e±*^-/2 /' " g{s)ds, u ^ ±oo. 



(4.12) 



Moreover, Karamata's theorem (Theorem 1.5.11(ii) in Bingham et al. (1987)) gives 



l/|a;| r-oo 

g{s)ds = s^ g{l/s)ds ~ \u!\^ g{l/\u}\)/h', uj — )• ±00, 

J\u)\ 

where we used the fact g{-) G 7^i^_i(0+) means ^(l/-) € 7^_,y+i(oo). 



Substituting (4.12) into (2.1) and recalling that T{v + 1) = uT{v), we obtain 



/y(|a;|) = ^\F{g)\\u:) ~ I^|a;rV(l/kl), ^ ^ ±00, 



which gives the desired result. 



□ 
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Remark 4.6. Condition {A2) can be replaced by a monotonicity condition on the derivative 



(^(M)(.) near the origin, so that the monotone density theorem (Bingham et al. (19871, Theo- 
rem 1.7.2.) can be apphed. D 



Example 4.7. [CARMA(p, q) process] 

The CARMA(p,q) process Y has spectral density (2.3), which clearly has the form 



/y(w) = |w| "L(|w|), U £ 



where a = 2{p — q) and lim^j-j^oo -^(I'^l) = '^^/(^t'")- Hence, by Theorem 4.2 ^c), the white noise 
variance in the Wold representation of Y satisfies as A J, 0, 



^2 ^2ri \2(p-q)-l 



(4.13) 



rem 



where C2(p_„) can be calculated from (4.4). However C2(p-q) "^^^ ^^^° ^^ calculated from (3.5) 
as C2(n-g\ = [(2(p - q) - l)!nf=f~ limA;o^(Ci)]~\ where ■q{^i) was defined in ( 2.16[ ). Theo- 



4.2 ^b) implies that the spectral density of A"^ P+'^/'^{\ — BY '^Y converges to that of a 



short memory stationary process. From Theorem 12.11 we get the more precise result that the 

1/2 



spectral density of C„, _ ^ 
noise with variance o"^. 



A'?-P+i/2(i _ B)P-<i nLi(l + v{ii)B)~^Y^ converges to that of white 

□ 



Example 4.8. [FICARMA(p, d, q) process, [Brockwell and Marquardt| ( |2005[ )] 
The fractionally integrated causal CARMA(p, d, q) process has spectral density 



frico) 



a' 1 

27r|a;|2'^ 



b{iuj) 
aiiuj) 



UJ £ Q.C 



(4.14) 



with a(-) and 6(-) as in (2.3) and < d < 0.5. Hence 



/y(w) = |L^r"L(|w|), U(^VL, 



where a = 2{p+d—q) and lim^^oo -^(I'^l) = '^ /(^''t). The spectral density (4.14) has a singularity 



at frequency which gives rise to the slowly decaying autocorrelation function associated with 



long memory. Applying Theorem |4.2[ c) as in Example 4.7, the white noise variance in the Wold 
representation of Y satisfies as A ^ 



2 Jlf-i A2(p+d-q)-l 



(4.15) 



where C2(p+d-q) can be calculated from (4.4). As A 4- 0, the asymptotic spectral density /a of y 



is given by (4.2) with a = 2{p + d — q) > 1 and is therefore not integrable for any A > 0. However 
Theorem 4.2 'b) implies that the spectral density of A'^~p^'^^^'^{1 — BY~^'^~'^Y^ converges to that 
of a short memory stationary process. □ 



Our next two examples are widely used in the modelling of turbulence. Kolmogorov's famous 



5/3 law (see Frisch (1996) Section 6.3.1, Pope (2000) Section 6.1.3) suggests a regularly varying 
spectral density model for turbulent flows. 
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Example 4.9. [Two turbulence models] 

Denote by U the mean flow velocity, with i the integral scale parameter and define i = i/U. 

(i) The von Karman ( 1948 ) spectrum models the isotropic energy spectrum. Its spectral density 

/ 9 \ 17/6 



is, for C and C£ positive, given by 



U) 



^w2 + ce/i 



LO G il.p. 



Moreover, fy € 7^-5/3, so it has a representation (4.1) and the conclusions of Theorem 
with a = 5/3. 



4.2 



hold 



(ii) The Kaimal spectrum for the longitudinal component of the energy spectrum is the current 
standard of the International Electrotechnical Commission; cf. lEC 61400-1 (1999). The spectral 
density is given by 



where v is the variance of Y. Moreover, fy G '^-5/31 so it has a representation (4.1) and the 
conclusions of Theorem 4.2 hold with a = 5/3. □ 



Example 4.10. [Gamma kernels and Whittle-Matern autocorrelations] 
The CMA process ( |1.1[ ) with gamma kernel, 

g(t) = t'^-ie-^*l(o,oo)(i), A > 0, u> 1/2, 



(4.17) 



has variance 

and autocorrelation function 

py{h) 



jy{0) = a'^{2X)^-^''r{2u - 1) 

r)3/2-U 



r(z.-l/2)''"^' "--v-vr-i^' (4.18) 

which is the Whittle-Matern autocorrelation function (see [Guttorp and Gneiting (2005)) with 



parameter z^ — 1/2, evaluated at Xh. The function -fCi^_i/2 in (|4.18 



of the second kind with index u — 1/2 (Abramowitz and Stegun 



is the modified Bessel function 



(1974), Section 9.6). 



Note that g G TZ,y-i{0+) and that it satisfies the assumptions of Proposition 4.5 From (2.1 ) with 
g{u}) = T{h'){X — iuj)^'^, we obtain the spectral density 






u 



-2u 



a^T^iy) 



uj G rtf.. 



2^((A/a;)2 + l)^ 

which belongs to 7^_2!/(oo) and slowly varying function L such that lim^j-j^oo L{ijj) = cr^r^(i/)/27r. 

Note that \i v = 5/6, then /y, like the von Karman spectral density of Example 4.9 (i), decays 
as w"^'^ for cj — )• 00, in accordance with Kolmogorov's 5/3 law. 



Theorem 4.2 gives the asymptotic form of the spectral density of the sequence {(1 — BYY^^ 
as A I 0, 



n& 



-"^ + (2vr)-^'^C {2u, l-^)+ {2n)-"'C {2u, 1 + ^ 



\UJ\ 



u) G i}fj. 



15 



The second-order structure function, 52(A) := E[(yA — ^)^]) plays an important role in the 
physics of turbulence. For the kernel ( 4.17[ ) with 1/2 < u < 3/2 its asymptotic behaviour as 
A 4 is given by 

52(A) = 27y(0)(l-py(A)), A>0, 



which, by the asymptotic behaviour as A ^ of K^_i/2{A) (see Abramowitz and Stegun (1974| 
Section 9.6), gives the asymptotic formulae. 



52(A) 

27y(0) 



2i-2'^f^^^^^(AA)2--i + 0{A% \l2<v< 3/2, 



i(AA)2|logA|+0(A3), 



I 4(1.-3/2) 



(AA)2 + 0(A 



2f-l^ 



V = 3/2, 
zv > 3/2, 



which can be found in Pope| ( [2000 ), Appendix G, and Barndorff-Nielsen et al. (2011). The first 
of these formulae can also be obtained as a special case of Corollary |4.4| with p = 1. □ 



5 Estimating the kernel function g 

Given observations of Y with A small, we estimate the kernel g by estimating the approximation 



g defined in (3.2), which converges by Theorem 3.2 for all CARMA{p,q) processes pointwise, 
as A J, 0, to ag, or to g if the driving Levy process is standardized so that Var(Li) = l. In all of 
the examples which follow we make this assumption since without it a and g are confounded. 
To estimate g it suffices to estimate the coefficients and white noise variance in the Wold 



representation (3.1) of Y , for which standard non-parametric methods are available. Being 



non-parametric they require no a priori knowledge of the order of the underlying CARMA 
process and moreover they can be applied to observations of any stationary causal CMA of the 
form (1.1). We now describe two useful algorithms for the estimation of g and apply them to 



simulated sampled data from a CARMA process, which has rational spectral density, and from 
a CMA process with non-rational spectral density. In Section [6] we estimate the kernel for a real 
data set without any parametric assumption. 



5.1 The algorithms 



The Wold coefficients ip^ and white noise variance o"^ appearing in the definition (3.2) of g 
are completely determined by the autocovariance function of the sequence Y and can be es- 
timated using the sample autocovariances of the observations of Y . These will be based on 
high-frequency observations of 1^ , i = 1, . . . ,n, with n very large and A(> 0) very small. The 
estimator of g (and of g) obtained by substituting the estimated Wold coefficients and white 



noise variance into (3.2) will be denoted g . 



The estimation of the coefficients ipj' for j G N, and the white noise variance a^ from the 
sample autocovariance function of Y is a time-domain version of the spectral factorization 



problem for which many algorithms have been developed (see Sayed and Kailath (2001) for a 
recent survey). 
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We shall apply and compare two algorithms for generating the required estimates. In both 
cases we fit high-order moving averages to Y . Note that g as defined by ( |3.2[ ) is piecewise 
constant on the intervals [jA, (j + 1)A) for j = 0, 1, 2, . . .. In Figures 2-5 and 6 c) the values of 
g are plotted at the mid-points (n -|- 1/2)A, for n = 0, 1, 2, . . .. 

5.1.1 Durbin-Levinson algorithm 

The first method determines the (high-order) causal AR(p) process whose autocovariances match 
the sample autocovariances of Y up to lag p by solving the Yule- Walker equations for the 
coefficients (f)i, . . . ,<j)p of the autoregressive polynomial (piz) = 1 — (f)iz — ... — (j)pZ^. This can 
be done efficiently using the Durbin-Levinson algorithm, which also yields an estimator of o"^ 



(Proposition 8.2.1, Brockwell and Davis (1991)). 



The coefficients of the Wold representation (3.1) of the fitted AR(p) process are used as 
estimates of the coefficients ipj in (4.1). The estimate of tjj^ is the coefficient of z^ in the power- 
series expansion of 1/0(2;). These coefficients are easily calculated recursively from the coefficients 



5.1.2 Innovations algorithm 

A more direct approach to fitting a high-order moving average process based on the sample 



autocovariance function is to use the innovations algorithm (Definition 8.3.1 of Brockwell and 



Davis (199ip). Although this method is computationally slower than the Durbin-Levinson al- 



gorithm, under the conditions of Theorem 8.3.1 of Brockwell and Davis (1991), the estimates 



of the Wold coefficients are asymptotically jointly normal with simple covariance matrix, and 



the estimator of the white noise variance o"^ is consistent. If the driving process L in (1.1) is 
Brownian motion (as in the following simulations) and if y is a CARMA process then Y and 
Y are Gaussian, the driving noise in the Wold representation of Y is iid Gaussian, and the 
conditions of Theorem 8.3.1 are satisfied. 

In practice it has been found that the Durbin-Levinson algorithm gives better results except 
when the fitted autoregressive polynomial has zeroes very close to the unit circle. 

5.2 Simulation examples 

We shall illustrate our estimation method by a small simulation study based on the gamma 



kernel of Example 4.10 with standard Brownian motion as driving noise process, A = 1, cr^ = 1 
and with the two values, ly = 1.05 and u = 2, the latter corresponding to the CAR(2) model 
with a{z) = (z — 1)^ and b{z) = 1. 



We applied the estimation method of Section 5.1 using both the Durbin-Levinson and the 
innovations algorithms. Since the models are known, we first applied both algorithms to the true, 
rather than the estimated autocovariances in order to eliminate the effect of sampling error. The 
estimated kernel functions are depicted in Figures 2-5 in the upper row. 

Then we simulated a sample path for each regime u = 1.05 and v = 2 based on n = 8 • 10^ 
time steps and took this as a proxy for the continuous time true sample paths. From this we 
considered two different time resolutions for the sampled process Y : A = 2~^ = 0.25 and 
A = 2~^ = 0.0625. We fitted the kernel function g{-) based on observations up to time T = 8, that 
is, we estimate g{{j + ^)A) for j = 0, . . . , iV = 8/2^2 = 32 and for j = 0, . . . , iV = 8/2-^ = 128, 
respectively. 
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Figure 2: Estimation of the gamnia kernel for u = 1.05 and A = 2 ^. 



In the case of the innovations algorithm, we used (for the known as well as for the estimated 
autocovariances) values of the discrete autocovariance functions up to 3A^. 

The results are depicted in Figures 2-5 bottom row, where the squares denote the estimates 
from the innovations algorithm, and the circles the ones from the Durbin-Levinson algorithm. 
For reference the true kernel function is plotted with a solid line. 

Comparing top and bottom rows of Figures 2-5 we find for the estimated autocovariance 
function an intrinsic finite-sample error, which influences the kernel estimation. We notice that 
in all cases considered, the Durbin-Levinson algorithm gives much better estimates. Further- 
more, as expected, the estimates for both algorithms improve with decreasing grid spacing. The 
Durbin-Levinson algorithm, unlike the innovations algorithm, provides estimates with are in good 
agreement with the original kernel function even for a coarse grid with A = 0.25. 



6 An application to real data: mean flow turbulent velocities 



We apply the algorithm of Section 5.1.1 to the Brookhaven turbulent wind-speed data, which 
consists of 20 x 10^ measurements taken at 5000Hz; i.e. 5000 data points per second. The series 
thus covers a total time interval of approximately 67 minutes. This dataset displays a rather high 
Reynolds number (about 17000), thus it can be regarded as a good representative of turbulent 
phenomena. A more detailed presentation of turbulence phenomena and an application of the 



CM A model ( 1.1 ) in the context of turbulence modelling is given in Ferrazzano and Kliippelberg 



(2012); moreover we refer to Drhuva (2000), Ferrazzano (2010) for a precise description of the 
data, and to Pope (2000), Frisch (1996) for a comprehensive review of turbulence theory. The 
non-parametric kernel estimate for the CMA model (1.1 ) is then based on A = 1/5000 = 2 x 10~^ 
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Figure 3: Estimation of the gamma kernel for u = 1.05 and A = 2 ^. 



seconds. A CMA model (1.1 ) with a gamma kernel as in Example 4.10 has been suggested as a 



parametric model in Barndorff-Nielsen and Schmiegel (2009) 



Figure Ig] a) shows the sample autocorrelation function up to 120 seconds, which appears to 
be exponentially decreasing. In general, the data are not significantly correlated after a lag of 
100 seconds. 

The estimated spectral density of Y is shown in Figure p^b), plotted against the frequency 
(p = ijj/2iT^ since uj in (2.1) have the units of an angular velocity. The estimates depicted with 
circles was estimated by Welch's method (Welch (1967)) with segments of 2^^ data points (circa 



14 minutes), windowed with a Hamming window and using an overlapping factor of 50%. This 
method allows a significant reduction of the variance of the estimate, sacrificing some resolution 
in frequency. In order to have a better resolution near to the frequency OHz, we estimated the 
spectral density for (^ < 10~^ Hz with the raw periodogram (Brockwell and Davis (1991 1, p. 322), 



which provides a better resolution in frequency at cost of a larger variance. The results are plotted 
in the leftmost part of Figure [6] b) with diamonds, and the two ranges of estimation are indicated 
by a vertical solid line. The spectral density is plotted on a log-log scale, so that a power-law 
relationship appears linear. The spectral density in a neighborhood of zero is essentially constant, 
and this is compatible with an exponentially decreasing autocorrelation function (as for instance 



for the gamma kernel function of Example 4.10). 

Moreover, for a frequency tp between 10~^ and 200Hz, log/y decreases linearly with log(/7 
with slope of approximately —5/3, in accordance with Kolmogorov's 5/3-law. As a reference, a 
solid line depicts the power-law ip~^'^. For tp larger than 200Hz, the spectral density deviates 
from Kolmogorov's 5/3-law, decaying with a steeper slope. We stress the point that a spectral 
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Figure 4: Estimation of the gamma kernel for z/ = 2 (CAR(2) process) and A = 2 . 



density decaying as prescribed by Kolmogorov's law in the neighborhood of oo would require a 
-1/6 near to the origin, according to Proposition 



4.5 



kernel behaving like t~ 

The estimated kernel function g {t) is plotted in Figure pic) in a log-linear scale, so that the 
behaviour of the kernel estimate for t near zero and for t to infinity is clearly visible. For large t 
the estimated g{t) decays rapidly and it oscillates slightly around zero for t > 100 seconds. For 
small t the estimated g{t) grows slowly up to t ~ 10~^, corresponding to Kolmogorov's 5/3-law. 
For smaller t it drops off to zero corresponding to the steeper decay of the spectral density at 
high frequencies. 

Figure 6 d) shows the spectral density computed directly from the estimated kernel func- 
tion g . Its close resemblance to the spectral density calculated by Welch's method provides 
justification for our estimator of g even when there is no underlying parametric model. 



7 Conclusions 

We studied the behaviour of the sequence of observations Y obtained when a CMA process of 
the form (1.1) is observed on a grid with spacing A as A ^ 0. 

In the particular case when y is a CARMA process we obtained a more refined asymptotic 
representation of the sampled process than that found by Brockwell et al. (2012) and used it to 
show the pointwise convergence as A J, of a sequence of functions defined in terms of the Wold 
representation of the sampled process to the kernel g. This suggested a non-parametric approach 
to the estimation of g based on estimation of the coefficients and white noise variance of the 
Wold representation of the sampled process. 
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Figure 5: Estimation of the gamma kernel for z/ = 2 (CAR(2) process) and A = 2 ^. 



For a larger class of CMA processes we found results analogous to those of Brockwell et al. 
(2012) and examined their implications for the local second-order properties of such processes, 
which include in particular fractionally integrated CARMA processes. 

Finally we applied the non-parametric procedure for estimating g to simulated and real data 
with positive results. 
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